library(gplots)

params <- matrix(nrow=1080,ncol=6)
dpl <- c(0,0.25,0.5,0.75,1)
scl <- c(0.001,0.01,0.1,0.5,1,2)
sml <- c(10^-5,10^-6,10^-7,10^-8)
nml <- c(10^-4,10^-5,10^-7)
sw <- 0

for (divprob in 1:5){
  for (selcoeff in 1:6){
    for (selmu in 1:4){
      for (neutmu in 1:3){
        for (replicates in 1:3){
          sw <- sw + 1
          params[sw,1] <- sw
          params[sw,2] <- dpl[divprob]
          params[sw,3] <- scl[selcoeff]
          params[sw,4] <- sml[selmu]
          params[sw,5] <- nml[neutmu]
          params[sw,6] <- replicates
        }
      }
    }	
  }
}

pdf("beneficial.mutations.frequencies.sm10-5.s1.pdf",width=8, height=6)
res <- matrix(data=NA,nrow=5,ncol=45)

column <- 0
for (r in 1:length(dpl)){
  for (nmu in 1:length(nml)){
    for (repl in 1:3){
      column <- column+1
      sw <- params[params[,3]==1&params[,4]==10^-5&params[,2]==dpl[r]&params[,5]==nml[nmu]&params[,6]==repl,1]
      f <- read.csv(paste("../sim",sw,"/sim",sw,".hit.log",sep=""),header=FALSE)
      f <- f[f[,1]==7300,]
      res[1:5,column] <- as.numeric(f[,3:7]) 
    }
  }
}
subres <- res[,1:9]
subres1 <- subres[,order(subres[1,]*1+subres[2,]*2+subres[3,]*4+subres[4,]*8+subres[5,]*16,decreasing=FALSE)]
subres <- res[,10:18]
subres2 <- subres[,order(subres[1,]*1+subres[2,]*2+subres[3,]*4+subres[4,]*8+subres[5,]*16,decreasing=FALSE)]
subres <- res[,19:27]
subres3 <- subres[,order(subres[1,]*1+subres[2,]*2+subres[3,]*4+subres[4,]*8+subres[5,]*16,decreasing=FALSE)]
subres <- res[,28:36]
subres4 <- subres[,order(subres[1,]*1+subres[2,]*2+subres[3,]*4+subres[4,]*8+subres[5,]*16,decreasing=FALSE)]
subres <- res[,37:45]
subres5 <- subres[,order(subres[1,]*1+subres[2,]*2+subres[3,]*4+subres[4,]*8+subres[5,]*16,decreasing=FALSE)]

nc <- t(t(c(NA,NA,NA,NA,NA)))
res <- cbind(subres1,nc,subres2,nc,subres3,nc,subres4,nc,subres5)

rownames(res) <- c("1 hit","2 hits","3 hits","4 hits", "5 hits")
scheme <- c("magenta","blue","green","orange","red")
res <- as.matrix(res)
mp <- barplot2(res, beside = FALSE,col = scheme[1:5],ylim = c(0,1), plot.grid = FALSE,axes=FALSE)
axis(2,at=c(0,.25,.5,.75,1),label=c(0,.25,.5,.75,1))
axis(1,at=c(c(5,15,25,35,45)*1.2-0.5),label=c("r=0","r=.25","r=.5","r=.75","r=1"),tick=FALSE)

#legend(24,1,legend=rownames(res),fill=gray(1-(0:4/4)))
mtext(side=2, text="Frequency of beneficial mutations",line=2)
mtext(side=1, text="Variation between individual runs",line=2)
dev.off()


